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COMPUTATION  OF  ELECTROMAGNETIC  SCATTERING  PARAMETERS 
FOR  LOGNORMAL  DISTRIBUTIONS  OF  MAGNETIC  SPHERES: 
THEORY  AND  ALGORITHMS 


1.  INTRODUCTION 

A  previous  report’  described  the  theoiy,  algorithmic  design,  and  testing  of  a 
subroutine  for  computing  the  electromagnetic  scattering  from  a  magnetic  sphere.  In  this 
report  the  methodology  for  using  this  subroutine  to  compute  the  mass  extinction,  scattering, 
absorption,  and  back^tter  coefficients  for  lognormally  distributed  particulate  ensembles  is 
developed. 


In  what  follows  the  relevant  parts  of  Mie  scattering  theonr  are  presented. 
Mass  extinction  coefficients,  and  the  lognormal  size  distribution  are  defined,  and  the  theory 
and  algorthims  for  integrating  Mie  scattering  parameters  over  size  distributions  are 
developed.  The  integrations  are  carried  out  in  terms  of  dimensionless  scarring,  and  size 
distribution  parameters,  which  are  simply  related  to  the  usual  mass  scattering  coefficients. 


2. 


LIGHT  SCATTERING  AND  CLOUD  MACROPHYSICS 


For  spheres,  the  efficiencies,  the  ratio  of  the  optical  cross  section  to  the 
geometric  cross  section,  for  extinction  (j2,),  scattering  (Q^),  absorption  (QJ,  and  backscatter 
(C»)  are:^ 


(2/1  + 1)  Re  (o„  -hbn  ) 


(1) 
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J.  +  Vi  (z).  H„  +  (z)  are  the  Bessel  functions  and  Hankel  functions  (second  kind)  of  half¬ 

integer  order.^  X  =tD/X  is  the  size  parameter  where,  D,  is  the  diameter  in  ;xm:  X,  the 

wavelength  in  /xm;  and  m  =  the  refractive  index  for  which  ]!,  the  permeability,  and 

r,  the  permittivity,  are  complex  quantities.  Primes  denote  derivatives  with  respect  to  the 
argument. 

For  ensembles  of  spherical  particles,  the  mass  extinction,  scattering, 
absorption,  or  backscatter  coefficients  (in  square  meters  per  gram)  at  wavelength  \  are 
defin^  as 


Q  =  S  /  E 
/  < 


(7) 


where  A,  is  the  optical  cross  section  for  extinction,  scattering,  absorption,  or  backscatter  of 
the  /th  particle  (in  square  meters),  is  the  mass  of  the  /th  particle  (in  grams),  and  the  sum 
extends  over  the  ensemble  of  particles  contained  in  the  optical  path.  For  a  continuous 
distribution  of  particles,  the  expression  for  a  becomes 


Q  *■  / Q{D)dm  ' 

a(P)  is  the  coefficient  of  extinction,  scattering,  absorption,  or  backscatter  for  a  particle 
of  size  D,  =  3Q  l2pD  where  Q  is  given  by  the  expressions  for  2,,  Q,,  Q„  and 
(^nations  1-4);  by  convention  an  additional  factor  of  4t  is  included  in  the  denominator  of 
this  expression  when  Q  —  Q^,.*  dm  is  the  mass-distribution  function  which  for  the  purposes 
of  this  work  is  chosen  to  be  the  lognormal  distribution  function: 


dm  =  (2;rln2  Cg  exp  {-V2[ln  (D  /D„)  /In  ag]-}d  {In  D) 


(9) 


3.  NUMERICAL  METHODS  -  ALGORITHM  DESIGN 

Evaluation  of  the  indicated  quadratures  over  lognormal  particle-size 
distributions  can  be  a  difficult  numerical  problem.  The  determination  of  the  2’-^  is  a 
computationally  intense  procedure  that  must  be  carried  out  over  a  fine  integration  mesh  to 
ensure  accurate  results;  this  is  especially  true  for  Q,  with  no  absorption  (Im(m)  =  0)  or  2* 
which  have  substantial  resonant  structure. 

Because  the  Q’s  are  evaluated  in  terms  of  the  dimensionless  size  parameter,  X, 
it  is  convenient  to  formulate  the  quadrature  problem  in  terms  of  dimensionless  quantities. 

The  a’s  can  be  put  in  a  dimensionless  form  as  follows; 
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Qj  =  Q  •  X  •  p 


(10) 


Where  p  is  the  density  in  g  cm'*.  Also,  converting  D,  D„  in  the  particle  size  distribution  to 
X,  X„  leaves  the  functional  form  of  the  distribution  unchanged. 

After  converting  these  quantities  to  dimensionless  form  the  integrations  are 
carried  out  by  establishing  minimum  (Xq)  and  maximum  {XD  size  parameters  that  are  the  end 
points  for  N  evenly  spaced  abscissae  on  a  logarithmic  grid  and  then  applying  an  appropriate 
numerical  integration  scheme  such  as  Simpson’s  rule.  The  X  grid  is  given  by 


Xo,  RXo. 


,R^-^Xo 


(11) 


where  Xo  =  XJc”  X,  =  =  R^'^  Xo,  and  /?  =  AlnX  =  n,  the  numter  of 

geometric  standard  deviations,  is  chosen  to  be  three  or  greater  so  that  regions  outside  the 
grid  contribute  a  numerically  negligible  amount  to  the  integral.  The  dimensionless 
coefficients  are  then  given  by  the  sum 


aj  =  {9n/S\n^<7gf^  ^  Xj  Qexp{-V2[{\nXj  -  InJ^^  )/lnajj  ]->  A(ln^)  ^^2) 


When  Q  =  Qb  then  the  expression  for  must  be  multiplied  by  an  additional  factor  of  llAr. 
The  Figure  shows  the  effect  of  varying  the  number  of  terms  in  the  sum  in  Equation  12  for 
the  dimensionless  extinction  coefficient  for  m  =  (3.0,  0.);  /i  =  (1.0,  0.). 


4.  RESULTS  AND  DISCUSSION 

The  listings  of  a  main  (driver)  program  (disint),  which  computes  either 
dimensioned  or  dimensionless  extinction,  scattering,  absorption,  and  back^tter  coefficients 
is  displayed  in  Appendix  A  along  with  the  listings  for  the  auxiliary  routines;  dxtcs,  magsph, 
zbjy,  and  simp.  The  Bessel  function  subroutine  package  by  Amos*  is  required  to  implement 
these  routines. 


Appendix  B  lists  the  input  to  disint  and  then  displays  the  resulting  output  for 
three  sample  calculations.  The  first  two  sample  problems  compute  the  dimensionless 
extinction  coefficients  for  a  mass  median  size  parameter  of  10,  a  refractive  index  of  (1.5, 
-0.1).  The  first  example  uses  a  very  narrow  distribution  (a,  =  1.01);  the  results  of  this 
calculation  can  be  inverted  by  the  relationship. 


<2  = 


2A'„ 


3r 
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(13) 
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and  the  results  compared  with  the  calculation  for  monodisperse  pa^cle  calculation  of 
Wiscombe*  (Table).  The  next  example  gives  the  results  for  a  particle  with  same  refractive 
index  and  a  geometric  standard  deviation  of  1.4.  The  final  example  is  a  sample  problem  that 
exercises  the  code  for  tiie  case  of  magnetic  spheres  with  /*  =  e  =  (2.24,  -.3),  which  are 
lognormally  distributed  with  an  MMD  of  10  ftm  and  a  ffg  of  1.5.  It  is  Imown  that  the 
backscatter  vanishes  for  spheres  with  these  material  properties;*"’  this  sample  problem 
illustrates  this  point. 


Table.  Efficiency  Factor  Comparisons 


Distribution 

Q. 

0. 

Q. 

Monodisperse 

2.4598 

1.2351 

1.2246 

Lognormal* 

2.4593 

1.2346 

1.2246 

\  =  1.01 


11 


Blank 


12 


LITERATURE  CITED 


1.  Milham,  Merrill  E.,  Electromagnetic  Scattering  by  Magnetic  Spheres: 
Theory  and  Algorithms.  ERDEC-TR-207,  U.S.  Army  Edgewood  Research,  Development  airf 
Engineering  Center,  MD,  October  1994,  UNCLASSIFIED  Report. 

2.  Van  de  Hulst.  H.C..  Light  Scattering  by  Small  Particles. 

Dover,  NY,  1981. 

3.  Abramowitz,  M.,  and  Stegun,  I.A.,  Handbook  of  Mathematical  Functions. 
GPO,  Washington,  DC,  1964. 

4.  Bohren,  C.F.,  and  Huffman,  D.R.,  Absorption  and  Scattering  of  Light  by 
Small  Particles.  Wiley,  New  York,  NY,  1983. 

5.  Amos,  D.E.,  "Algorithm  644:  A  Portable  Package  for  Bessel  Functions  of 
a  Complex  Argument  and  Nonnegative  Order,"  ACM  Trans  Math  Software  Vol  12, 

pp  265-273,  1986. 

6.  Wiscombe,  W.J.,  Mie  Scattering  Calculations:  Advances  in  Technique  and 
Fast.  Vector-Speed  Computer  Codes.  NCAR/TN-140+STR,  National  Center  for 
Atmospheric  Research,  Boulder,  CO,  1979. 

7.  Kerker,  M.,  Wang,  D.-S.,  and  Giles,  C.L.,  "Electromagnetic  Scattering  by 
Magnetic  Spheres,"  J.  Opt.  Soc.  Am.  Vol.  73,  pp  765-767  (1983). 


13 


Blank 


14 


APPENDIX  A:  PROGRAM  LISTINGS 


program  disint 

c************************************************************************* 

c 

c  Program  disint  computes  extinction,  scattering,  absorption,  and 
c  backscatter  coefficients  for  lognormal  distributions  of  homogeneous 
c  magnetic  or  nonmagnetic  spheres, 
c 

c  Merrill  Milheun  »>  version  1.0  «<  JANUARY  1994 

c 

c  inputs : 

c  ALL  INPUT  IS  IN  LIST  DIRECTED  FORMAT 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


line  #1 

mmd  =  mass  median  diameter  (micron)  (real*8) 
sigmag  «  geometric  standard  deviation  (real*8) 

rho  «  density  of  the  particulate  material 
if  rho  .le.  zero  dimensionless 
extinction,  scattering,  absorption 
and  backscatter  coefficients  are 
computed,  (g  cm-3)  (real*8) 

ndfp  B  number  of  logarithmically  equally  spaced 
abscissa  values  over  which  the  various 
scattering  coefficients  are  to  be 
integrated  (integer) 

nsig  =  number  of  sigmas  on  each  side  of  the  mmd  over 
which  the  scattering  coefficients  are  to  be 
integrated.  No  correction  is  made  for  excluded 
distribution  tails.  (real*8) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


line  #2 

nwl  -  number  of  wavelengths  and  material 

properties  to  be  read  (integer) 

flag  =  'r'  for  refractive  index  data  or 

'p'  for  permittivity  data  (character*!) 

wl  -  wavelength,  up  to  10  values  (micron)  (real*8) 
mp  =  complex  index  or  permittivity  values 

to  be  read,  up  to  10  values  (complex*16) 

line  #3 

muu  complex  permeability  values  to  be  read 

nwl  values  are  required  (complex*16) 


c 

c  output- 

c  line  1:  mmx  or  mmd  (micron)  sigmag  density  (g  cm-3)  ndfp  nsig 
c  line  2:  wavelength  (micron)  ext.  coef.  (m  sq  /  g)  scat.  coef.  (m  sq  /  g) 
c  line  3:  abs.  coef.  (m  sq  /  g)  bsca.  coef.  (m  sq/g/sr) 

c  line  4:  permittivity  permeability 
c  line  5:  refractive  index 


c 

c  If  dimensionless  quantities  are  computed,  the  output  is  so 

c  labeled. 


c 

c  subroutines  used: 
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c 

c  dxtcs  -  returns  integrated  extinction  coefficients  for  a 

c  homogeneous  magnetic  sphere  or  homogeneous 

c  nonmagnetic  sphere  if  the  permeability  >  (1/0) 

c 

Q*  W  iNriAr  ***  ******  ♦***«*♦****★  ***  4nNr  ★«***  ******  ♦«**«**  ****  «nlnllr 

C 

implicit  none 
real*8  pi 

parameter  (pi*3 . 14159265358979d0) 
real*8  mmx/sigmag/nsig,dxe/dxs/dxa/dxb 
real *8  mmd,xC/rho,wl(10) /Xe/xs,xa,xb 
Integer  ndf p , i / nwl 
complex*16  muu(lO) /mU/mp(10) /eps 
real*8  one, zero 
parameter  (one-1 .dO, zero^O.dO) 
character*!  flag 
logical  dflag 
c 


write (*/*)  'read(*,*)  mmd, sigmag, rho,ndfp, nsig' 
read(*,*)  mmd, sigmag, rho, ndf p, nsig 

write(*,*)  'read(*,*)  nwl , flag, (wl (i) ,mp(i) , i=l , nwl ) ' 
read(*, *)  nwl, flag, (wl(i) ,mp(i) ,i=l,nwl) 
write(*,*)  'read(*,*)  (muu(i) ,i=l,nwl) ' 
read(*,*)  (muu(i) , i=l,nwl) 
write (*/ *) 


c 

* 


* 


df lag=rho. le . zero 
if (dflag)  then 
rho=one 
write (*, 1050) 

write (*, 1100)  mmd, sigmag, rho, ndfp, nsig 
else 

write (* , 1000)  mmd, sigmag, rho, ndfp, nsig 
end  if 
write(*, *) 

do  100  i=l,nwl 


if(dflag)  then 

mmx=mmd 

wl  ( i ) =one 

mmx=pi *mmd/wl ( i ) 
end  if 


else 


if (flag.eq. 'p' )  then 
eps=mp(i) 

else  if (flag.eq. 'r' )  then 
eps=(mp(i) ) **2 

else 

stop  'material  properties  input  error' 
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end  if 


mu^uu  ( i ) 

call  dxtcs (mmx , sigmag , ndf p , ns ig , eps , mu , dxe , dxs , dxa , dxb ) 

xc=l . do /wl ( i ) /rho 

xe=xc*dxe 

xs=xc*dxs 

xa=xc*dxa 

xb=xc*dxb 

write(*,1200)  wl(i) ,xe,xs,xa,xb 
write  (*,1300)  eps, mu 
if (flag. eg. 'r' )  write(*,1400)  mp(i) 
write (*,*) 

100  continue 

200  stop  'done  disint' 

1000  format ('mmd  *  ',f5.2,'  sigmag  =  ',f5.2,'  rho  =  ',f5.2, 

&  '  ndfp  =  ',i5,'  nsig  =  ',f5.2) 

1050  format ( 'dimensionless  extinction  coefficients  will  be  computed'/) 
1100  format ('mmx  =  ',f5.2,'  sigmag  =  ',f5.2,'  rho  *  ',f5.2, 

&  '  ndfp  *  '/is,'  nsig  =  ',f5.2) 

1200  format('wl  =  ',f8.4,'  ext  =  ',el0.5,'  scat  =  ' ,el0.5, /13x, 

&  '  abs  —  ',el0.5,'  bsca  =  ',ei0.5) 

1300  format ('eps  =  ',2f6.3,'  mu  =  ',2f6.3) 

1400  format (' refractive  index  =  ',2f6.3) 

end 
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, dxe , dxs , dxa , dxb ) 


subroutine  dxtcs (mmx , sigmag, ndfp, nsig, eps ,mu 


c»>  subroutine  dxtcs  computes  dimensionless  extinction,  scattering,  <« 
c»>  absorption,  and  backscatter  coefficients  for  particulate  ensembles  «< 
c»>  of  spheres  with  lognormal  particle  size  distributions  «< 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

* 


Merrill  Mllham 


»>  version  1.0  «< 


JANUARY  1994 


Inputs : 


mmx  mass  median  size  parameter  (real*8) 

geometric  standard  deviation  (real*8) 

nu^er  of  equally  spaced  (logarithmically) 
points  at  which  integrand  is  to  be 

evaluated  (integer) 

number  of  geometric  standard  deviations  on 
both  sides  of  the  mmx  over  which  the  integration 
extends  (real*8) 

complex  permittivity  of  particle  material  (complex*16) 
complex  permeability  of  particle  material  (complex*16) 


sigmag 

ndfp 


nsig 


eps 

mu 


outputs : 


dxe  “ 
dxs  - 
dxa  B 
dxb  = 

subroutines  used: 


dimensionless  extinction  coefficient  (real*8) 
dimensionless  scattering  coefficient  (real*8) 
dimensionless  absorption  coefficient  (real*8) 
dimensionless  backscatter  coefficient  (real*8) 

simp(real*8  function  )  -  does  Simpson's  rule 
integration 


implicit  none 


real*8  mmx, sigmag, nsig 
integer  ndfp 
complex*16  eps, mu 
real*8  dxe, dxs, dxa, dxb 

integer  mpsd 
parameter  (mpsd=3001) 

real*8  tsume (mpsd) ,tsums (mpsd) ,tsumb(mpsd) ,dx, dc 

real*8  Insg, e22x, r, xo, zo, z 

real*8  theta, qext, qsca, qbac, g 

real *8  simp 

complex*16  sl,s2 

integer  i 


real*8  pi 
real*8  dco,sqrt2 

parameter  (pi=3 . 14159265358979d0) 
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parameter  (dco=1.5d0*1.25331413731550d0,sqrt2*1.41421356237309d0) 
realms  zero,one, two, fourpi 

parameter  (zero=0.d0,one*l.d0,two=2 .dO, fourpi*4.d0*pi) 

Insgsdlog(sigmag) 

dc>dco/lnsg 

r»sigmag** (two*nsig/dble (ndfp-1 ) ) 

xo«mmx/(sigmag**nsig)/r 

zoBone/sqrt2/lnsg 

90  if (ndfp.gt .mpsd)  stop  'dxtcst  arrays  too  small  $90' 

theta-zero 
do  100  iBl,ndfp 
xo*xo*r 

zszo*dlog ( xo/mmx ) 
ez2x*dexp(-z*z)/xo 

call  magsph ( xo , eps , mu , 0 , theta , qext , qsca , qbac , g , s 1 , s2 ) 

tsume (i)sqext*ez2x 
tsums (i)Bqsca*ez2x 
tsumb(i)«qbac*ez2x 

100  continue 

dxsdlog(r) 

dxe-dc*simp ( tsume , ndfp, dx ) 
dxsBdc*simp ( tsums , ndfp, dx) 
dxa-ddim ( dxe , dxs ) 
dxb»dc / f ourpi *  s imp ( t sumb , ndf p , dx ) 

return 

end 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  magsph ( x , eps , mu , numang , theta , 

&  gext,qsca,qbac,g,sl,s2) 


Subroutine  magsph  computes  the  scattering  cross  sections  and  angular 
scattering  from  a  magnetic  sphere.  If  the  number  of  scattering  angles 
is  set  to  zero,  only  the  cross  sections  (efficiencies)  are  returned. 

Merrill  Milheun  »>  version  2.0  «<  SEPT  1993 


Inputs: 

X  «  size  parameter  of  the  sphere 
eps  «  complex  permittivity:  epsr  -i*epsi 
mu  «  complex  permeability:  mur  -  i*mui 
numang  *  number  of  scattering  angles 
between  0  &  90  deg. 
theta  =  scattering  angles  in  degrees 

theta (i)  are  entered  between  0  &  90  deg. 
theta  must  increase  monotonically .  Results  for 
supplementary  angles  (180  deg.  -  theta(i))  are 
also  returned. 


(real*8) 

(complex*16) 

(compl6X*16) 

(integer) 

(real*8) 


Outputs : 

qext  =  extinction  efficiency 
qsca  «=  scattering  efficiency 
qbac  =  backscatter  efficiency 
g  *  asymmetry  factor 
si  =  scattered  amplitude 
s2  «  scattered  amplitude 


(real*8) 

(real*8) 

(real*8) 

(real*8) 

(complex*16) 

(complex*16) 


Subroutines  used : 


zbjy  -  returns  one-half  integer  order  J  &  Y  Bessel  functions 
References : 


M.  Kerker,  D.-S.  Wang,  and  C.  L.  Giles,  "Electromagnetic  scattering 
by  magnetic  spheres,"  J.  Opt.  Soc.  Am.,  73,  765-767  (1983). 

D.  E.  Amos,  "Algorithm  644: A  portable  package  for  Bessel  functions  of 
a  complex  argument  and  nonnegative  order,"  ACM  Trans,  on  Math. 
Software,  12,  265-273  (1986). 

M.  Abramowitz  and  I.  A.  Stegun,  "Handbook  of  Mathematical  Functions," 
NBS  Applied  Math.  Series  55,  US  Dept,  of  Commerce,  Washington,  DC 
(1955). 


W.  J.  Wiscombe, "Hie  Scattering  Calculations:  Advances  in  Technique  and 
Fast,  Vector-Speed  Computer  Codes,"  NCAR  Tech.  Note,  NCAR/TN-140+STR 
(1979) 
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c 

c 

c 

it 

C 

if 


it 

it 


C 

* 


c 


icificitieicicititirirititititiritititiciciticititititifieiticiticititiriricieiriririticieitieieicifiticicirititititicifitititirititititititiritit 

implicit  none 
real*8  x 

complex*16  eps,mu 
integer  numang 
real*8  theta(l) 

real*8  gext,qsca,qbac,g 
complex*16  sl(l),s2(l) 

integer  al,nangl,nangl2 
real*8  third 

parameter  ( third=l . dO/3 . dO , al=5100 , nangl=255 , nangl2= (nangl+1 ) /2 ) 
complex*16  sp(nangl2) ,sm(nangl2) ,sps(nangl2) ,sms(nangl2) 

complex*16  m,mcl,mxi,s,t,u,v,an,bn,xp 
real*8  xi,dn,dnn,rn,tnpl,thetan 

real*8  xmu(nangl) ,pi(nangl) ,pil(nangl) ,tau(nangl) 
real*8  bjr (al) ,byr (al) 

real*8  cjr(al) ,cji(al) ,cyr(al) ,cyi(al) 

real*8  bjn,byn,bjl 

real*8  sc,ca,ti,t2,t3,t4 

complex*16  cjn,cyn,cjl,cyl,h2n,h21 

complex*16  anl , bnl , bs , anp, bnp, abp, abm, anpm, bnpm 

complex*16  ztl,zt2 

integer  kstop,k,mm,n, j , j2, j3 

real*8  cpi, zero, one, two, rad, half , fnu 
complex*16  cdblei,cdblel,cdbleO 

parameter  (cpi=3.1415926535897932384d0,zero=0.d0,one=l.d0) 
parameter  (two=2 .dO, rad=cpi/ 180 .dO, half =0 . 5d0, fnu=half) 
parameter  (cdblei=(0.d0, l.dO) ,cdblel=(l.d0,0.d0) ) 
parameter  (cdble0=(0.d0,0.d0) ) 

kstop=idint (x+4 .d0*x**third+4.d0) 

if (kstop.le.al)  then 
continue 

else 

print*, 'magsph  arrays  too  small:  kstop  =',kstop,'  al  “',al 
stop 
end  if 


if (numang. eq.O)  then 
si (l)=cdble0 
s2 (l)=cdble0 

else 


Appendix  A 


21 


if (numang.le.nangl2)  then 
continue 

else 

print*, numang, 'scattering  angles  input:  only ' ,nangl2, '  allowed' 
stop 
endif 

do  100  n=l,numang 
thetan^dabs (theta(n) ) 
theta (n)«thetan 

if (thetan.le.90.d0)  then 
continue 

else 

print*, 'theta( ' ,n, ')=' ,thetan, 'scattering  angles  must  be  <  90  deg' 
stop 
end  if 

thetan=rad*thetan 
xmu (n)=dcos (thetan) 

sp(n)=cdbleO 
sm(n)=cdbleO 
sps (n)=cdble0 
sms (n)=cdble0 
pi(n)‘=half 
pil(n)  =zero 

100  continue 

end  if 

m=2sqrt (mu*eps ) 
mcl=m/mu 

xi=one/x 

mxi-xi/m 


call  zbjy (x,m,kstop, fnu,bjr,byr,cjr,cji) 
bjl=bjr(l) 

c j l=dcmplx (c j  r ( 1 ) , c j i ( 1 ) ) 
cyl=dcmplx(cyr(l) ,cyi(l) ) 
h21=dcmplx(bjr(l) ,-byr(l) ) 


qexts=zero 

qsca=:zero 

bs^cdbleO 

g=zero 

anl=cdble0 

bnl=cdble0 
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do  300  k-2,kstop 


tnpl=tnpl+two 

tl=dn-rn 

ca=one+rn 

sc=rn 

dnn=dn+one 

rn^one/dnn 

sc®sc+rn 

bj n=bjr (k) 
byn=byr(k) 

cjn=dcmplx(cjr (k) ,cji (k) ) 
cyn=dcmplx(cyr (k) ,cyi (k) ) 
h2n»dcmplx (b j  n , -byn ) 


xp=dn*mxi 

s=cjl-xp*cjn 

u=s*h2n 

s=s*bjn 

xp-dn*dcmplx(xi, zero) 
t=cjn* (bjl-xp*bjn) 
v=cjn* (h21-xp*h2n) 

an®(s-mcl*t) /(u-mcl*v) 
bn=(nicl*s-t )  /  (mcl*u-v) 
abp=an+bn 
abm=an-bn 

ztl=dconjg(an) 
zt2=dconjg(bn) 
qext=qext+tnpl*dble(abp) 
qsca=qsca+tnpl* (an*ztl+bn*zt2 ) 
bs-bs- (dn+half  )  *inm*abm 

g=g+tl*dble(anl*ztl+bnl*zt2)+sc*dble (an*zt2) 

if (numang.eq.O)  then 
continue 

else 

anp=sc*abp 

bnp=sc*abm 

anpm=nmi*anp 

bnpm=inm*bnp 
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do  375  j-l,numang 

tl=xmu(j)*pi(j) 

t4=tl-pil(j) 

tau(  j  )*=dn*t4-pil  ( j  ) 

t2=pi(j)+tau(j) 

t3=pi(j)-tau(j) 

sp(j)=sp(j)+anp*t2 
sms ( j )=sms ( j )+bnpm*t2 
sm( j )=sm( j )+bnp*t3 
sps(j  )»=sps(j  )+anpm*t3 

Pil(j)=pi(j) 
pi ( j )=tl+ca*t4 

375  continue 
end  if 

dn=dnn 

nun^-mm 

anl^an 

bnl^^bn 


bjl=bjn 

cjl=cjn 

cyl^cyn 

h21=h2n 

300  continue 

if (numang.eq.O)  then 
continue 

else 

j2=2*numang 

do  500  j=l,numang 

j3=j2-j 

sl(j)=sp(j)+sm(j) 
s2(j)=sp(j)-sm(j) 
sl(j3)=sps(j)+sms(j) 
s2(j3)=sps(j )-sms(j ) 
500  continue 

end  if 


xi=two*xi*xi 

qext=xi*qext 

qsca=xi*qsca 

xi*two*xi 

qbac=xi *bs  *dcon j  g (bs ) 
g=xi/qsca*g 

return 
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subroutine  zbjy  (x,in,nstop,  fnu, 

&  bjr,byr,cjr,cji) 

c 

c  ************************************************************************* 

c 

c  Subroutine  zbjy  gets  J  &  Y  Bessel  functions  for  use  in 
c  sphere  (fnu=0.50)  or  cylinder  (fnu=O.OdO)  scattering  calculations, 
c  Scaled  or  nonscaled  functions  for  argument  z  *  m^x  are  returned  depending 
c  on  the  magnitude  of  the  product  of  the  size  parameter  and  the  complex 
c  refractive  index  (zabs(z)).  Nonscaled  functions  are  returned  if  the 
c  imaginary  part  of  the  refractive  index  is  zero.  Nonscaled  functions  are 
c  returned  for  argument  x. 
c 

c  Merrill  Milham  »>  version:  2.0  «<  JANUARY  1994 

c 

c  inputs: 

c  X  =  the  size  parameter  of  the  cylinder. (real*8) 

c  m  =  the  complex  refractive  index,  n  -  ik. (complex*16) 

c  nstop  =  the  highest  order  of  the  Bessel  functions . (integer) 

c  fnu  «=  O.SdO  for  sphere  calculations  or 

c  O.OdO  for  cylinder  calculations.  (real*8) 

c 

c  outputs : 
c 

c  bjr  *  real  part  of  J(x)  Bessel  functions . (array :  real*8) 

c  byr  “  real  part  of  Y(x)  Bessel  functions . (array :  real*8) 

c  cjr  *  real  part  of  J(m*x)  Bessel  functions . (array :  real*8) 
c  cji  *  imag.  part  of  J(m*x)  Bessel  functions . (array :  real*8) 
c 

c  subroutines  used: 
c 

c  zbesj  -  returns  J  Bessel  functions 

c  zbesy  -  returns  Y  Bessel  functions 

c 

c  Reference:  D.  E.  Amos,  "Algorthim  644:  A  Portable  Package  for  Bessel 
c  Functions  of  a  Complex  Argument  and  Nonnegative  Order, " 

c  ACM  Transcations  on  Mathematical  Software,  12,265-273(1986). 

c 

c**************************************************************************** 

c 

implicit  none 

automatic  cwrkr,cwrki,bji,byi 

«r 

integer  al 
real*8  zero, xll, two 

parameter  (al=5100, zero=0 . OdO, xll=l .d-1 , two=2 .dO) 

* 

real*8  bjr(l) ,byr(l) 
real*8  x, fnu, cjr (1 ) ,cj i (1) 
real*8  zr,zi,dlmach 
real*8  cwrkr (al) ,cwrki (al) 
real*8  rlm5,elim,aa,alim 
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coinplex*16  m,z 

integer  kode,ierr,nz,nstop,n 
integer  kl, iimach,k2,k 
logical  zfiag,eflgl,eflg2 

kode-1 

call  zbesj (x,zero, fnu,kode,nstop,bjr,cji,nz,ierr) 

eflgl^ierr.eg.O 
eflg2‘=nz.eq.O 
if  (eflgl.and.eflg2)  then 
continue 

else  if  ( .not .eflg2 .and.eflgl)  then 
nstop=nstop-nz 

print*, 'zbesj  error:  ierr  **',ierr,'nz  «',nz 
else 

print*, ' inputs  =' , x, zero, fnu,kode,nstop 
print*, 'zbesj  called  from  subroutine  zbjy' 

end  if 

call  zbesy (X, zero, f nu, kode , nstop, byr , c j i , nz , cwrkr , cwrki , ierr) 

eflgl=ierr.eq.O 
eflg2=nz .eq.O 
if  (eflgl.and.eflg2)  then 
.  continue 

else  if  ( .not. eflg2. and.eflgl)  then 
nstop=nstop-nz 

print* ,' zbesy  error:  ierr  =',ierr,'nz  »'',nz 

continue 

else 

print*, ' zbesy  error:  ierr  =',ierr,'nz  =',nz 
print*, ' inputs  =' , x, zero, fnu, kode, nstop 
print* ,' zbesy  called  from  subroutine  zbjy' 

end  if 

z=m*x 

zr=dble(z) 

zi=dimag(z) 

if (zi.ne.zero)  then 

kl  =  ilmach(15) 

k2  »  ilmach(16) 

rlmS  s  dlmach(5) 

k  «  min0(iabs(kl),iabs(k2)) 

elim  *  2.303d0*(dble(float(k))*rlm5-3.0d0) 

kl  ®  ilmach(14)  -  1 

aa  rlm5*dble(float  (kl) ) 

aa  —  aa*2.303d0 

alim  =  elim  +  dmaxl(-aa,-41.45d0) 
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if (dabs (zi ) .gt . alim)  kode=2 


100 


200 


continue 
end  if 


else 


call  zbesj (zr,zi,fnu,kode,nstop,cjr,cji,nz,ierr) 

eflgl=ierr .eq.O 
eflg2=nz .eq.O 
if  (eflgl .and.eflg2)  then 
continue 

else  if  (. not. eflg2. and. eflgl)  then 
nstop=nstop-nz 

print*, 'zbesj  error;  ierr  *',ierr,'nz  ■ 

continue 

else 

print*, ' zbesj  error;  ierr  =',ierr, 'nz  *',nz 
print*, 'inputs  =' , zr , zi , fnu,kode,nstop 
print *,' zbesj  called  from  subroutine  zbjy' 

end  if 


zflag=dabs(zi) .lt.two*dlmach(l) .and.x.lt.xll 
if  (.not.zflag)  then 
continue 

else 

do  100  n=l,nstop 
cj i (n)=zero 
continue 
end  if 

zflag=dabs(zr) .It. two*dlmach(l) .and.x.lt.xll 
if  (.not.zflag)  then 
continue 

else 

cji (l)=zero 
do  200  n=2,nstop 

if (mod(n,2) .eq.O)  then 
cjr (n)=zero 

else 

cji(n)=zero 
end  if 
continue 
end  if 

return 

end 
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real*8  function  simp  (y,n,dx) 

c***************************************************************** 

function  simp  does  simpson's  rule  integration 
inputs : 

y  *  array  containing  ordinate  values 
n  number  of  ordinate  values 
dx  abscissa  increment  (constant)  between 
ordinate  values 

output : 

simp  estimated  value  of  the  integral 


(real*8) 

(integer) 

(real*8) 

(real*8) 


c 
c 
c 
c 
c 
c 
c 
c 
c 

Q*********** ***********************************  ******************* 

c 

implicit  none 
real*8  si, s2,s4, third 
real*8  dx,y(l) 
integer  n,ne,i 
parameter  (third=l.d0/3.d0) 
c 

if (n. It. 5)  then 
write (6, 1000) 

1000  format ('  error  function  :  simp') 

write(6,'(''  number  of  integration  points  must  be  .ge.  5  n  « 
+  ,il0)')  n 

stop 
c 

else  if(mod(n,2)  .ne.l)  then 
write (6, 1000) 

write(6,'(''  number  of  integration  points  must  be  odd.  n  ■  ' 
+  ' )  n 

stop 
c 

else 

continue 
end  if 
c 

ne=n-3 

sl=y(l)+y(n) 

s2=0.d0 

s4=y(n-l) 

c 

do  100  i=2,ne,2 

s4=s4+y(i) 

s2*s2+y(i+l) 

100  continue 


c 

c 


s4-4.d0*s4 
s2»2 .d0*s2 

simp=third*dx* (sl+s2+s4) 

return 
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rilO) 


Blank 
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APPENDIX  B:  SAMPLE  OUTPUT 


INPUT 


10,1.01,-1,151,5 
1 , ' v' , 1 , (1 . 5 , — • 1) 
(IrO) 


OUTPUT 


read(*,*)  amd,sigmag,rbo,ndfp,nsig 
read(*,  *)  nwl,flag,  (wl(i)  ,mp(i) ,i’^l,nwl) 
read(*,*)  (muu(i) ,i^l,awl) 

dimensionless  extinction  coefficients  will  be  computed 

mmx  »»  10,00  sigmag  1.01  rbo  **  1.00  ndfp  »  151  nsig  * 

wl  =  1.0000  ext  =  .11589E+01  scat  «  .58179E+00 

abs  «  .57709E+00  bsca  «  .33692E-02 
eps  «  2.240-0.300  mu  »  1.000  0.000 

refractive  index  =  1.500-0.100 


INPUT 


10,1.4,-1,151,5 
1 ,  '  r ' ,  1 ,  (1 . 5 ,  — .  1 ) 

(IfO) 

OUTPUT 


read(*,*)  mmd, sigmag, rbo, ndfp, nsig 
read(*,*)  nwl,flag,(wl(i),mp(i),i=l,nwl) 
read(*,*)  (muu(i)  ,  i=l,nwl) 

dimensionless  extinction  coefficients  will  be  computed 

mmx  -  10.00  sigmag  =  1.40  rbo  1.00  ndfp  —  151  nsig  = 

wl  «  1.0000  ext  =  .12361E-h01  scat  *  .63179E+00 

abs  «  .60428E+00  bsca  *  .35034E-02 
eps  s  2.240-0.300  mu  1.000  0.000 
refractive  index  *»  1.500-0.100 


INPUT 


10,1.4,-1,151,5 
1 ,  'p ' , 1 , (2 . 24 ,  — . 3) 
(2. 24, -.3) 


5.00 


5.00 


OUTPUT 
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read (*,*)  mad,  sigmag,  rho,  ndfp,  nsig 
read(*,  *)  nwl,flag,  (wl  (i)  ,mp(i)  ,  i=J, nwl) 
read(*,  *)  (muu(i)  ,i’=l,nwl) 

dimens ioDless  extinction  coefficients  will  be  computed 

mmx  «  10.00  sigmag  «  1.40  rbo  «  1.00  ndfp  «  151  nsig  « 

wl  «  1.0000  ext  «  .12292E+01  scat  *  .57928E-h00 

abs  *  .64991E+00  bsca  *  .42778E-32 
eps  *  2.240-0.300  mu  *  2.240-0.300 


5.00 
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